Data and Libraries

Libraries

  library(phyloseq) #if not installed: install.packages('devtools')
  library(devtools) #if not installed: install_github("microbiome/microbiome") # Install the package
## Loading required package: usethis
  library(microbiome)
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.5.2
## 
## microbiome R package (microbiome.github.com)
##     
## 
## 
##  Copyright (C) 2011-2022 Leo Lahti, 
##     Sudarshan Shetty et al. <microbiome.github.io>
## 
## Attaching package: 'microbiome'
## The following object is masked from 'package:ggplot2':
## 
##     alpha
## The following object is masked from 'package:base':
## 
##     transform
  library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
  library(ggplot2)
  library(cowplot)
  library(RColorBrewer)
  library(remotes) #if not installed: remotes::install_github("jfq3/ggordiplots")
## 
## Attaching package: 'remotes'
## The following objects are masked from 'package:devtools':
## 
##     dev_package_deps, install_bioc, install_bitbucket, install_cran,
##     install_deps, install_dev, install_git, install_github,
##     install_gitlab, install_local, install_svn, install_url,
##     install_version, update_packages
  library(ggordiplots)
## Loading required package: vegan
## Loading required package: permute
## 
## Attaching package: 'permute'
## The following object is masked from 'package:devtools':
## 
##     check
## 
## Attaching package: 'vegan'
## The following object is masked from 'package:microbiome':
## 
##     diversity
## Loading required package: glue
  library(vegan)
  library(ALDEx2) #to install: BiocManager::install("ALDEx2")
## Loading required package: zCompositions
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## Loading required package: truncnorm
## Loading required package: survival
## Loading required package: lattice
## Loading required package: latticeExtra
## 
## Attaching package: 'latticeExtra'
## The following object is masked from 'package:ggplot2':
## 
##     layer
  library(decontam)
  library(brms)
## Loading required package: Rcpp
## Loading 'brms' package (version 2.23.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
## 
## Attaching package: 'brms'
## The following object is masked from 'package:survival':
## 
##     kidney
## The following object is masked from 'package:phyloseq':
## 
##     nsamples
## The following object is masked from 'package:stats':
## 
##     ar
  library(MoMAColors)
  library(ggrepel)
  library(ggvenn)
  library(gllvm)
## Loading required package: TMB
## 
## Attaching package: 'gllvm'
## The following object is masked from 'package:vegan':
## 
##     ordiplot
## The following object is masked from 'package:stats':
## 
##     simulate
  library(microbiomeutilities) #devtools::install_github("microsud/microbiomeutilities")
## Warning: replacing previous import 'ggplot2::alpha' by 'microbiome::alpha' when
## loading 'microbiomeutilities'
## 
## Attaching package: 'microbiomeutilities'
## The following object is masked from 'package:microbiome':
## 
##     add_refseq
  library(ggvegan)

Global Plotting Options

Quick access options for output graphics to make legends and axis labels larger / more legible

### Site Colours for Visit One From Paired Palette
  sitecols<-brewer.pal(8,"Paired")[c(1,3,5,7)]

Data

         ####### READ IN FILES
          
            #Taxonomy
              tax1<-readRDS('Newt 2021 Taxonomy.rds')
              #rownames(tax1)<-paste0("SV",seq(nrow(tax1)))
          
            #Sample Metadata  
              sample1<-readRDS('Newt 2021 SampleData.rds')
              head(sample1)
##   Sample
## 1   DP-1
## 2  DP-10
## 3  DP-11
## 4  DP-12
## 5  DP-13
## 6  DP-14
              colnames(sample1)[1]<-"swab"
              sample1$visit<-ifelse(grepl("C$",sample1$swab),2,1)
              sample1$swab_time<-with(sample1,paste(gsub("C$","",swab),visit,sep="_"))
              write.csv(sample1,'Newt Sample Data.csv')
              
            ## Additional Metadata
              meta1<-read.csv('Newt Metadata Full Feb17 R Input.csv')
              head(meta1)
##   visit airtemp watertemp   ph sex svl infection swab
## 1     1    12.3      13.1 5.11   1 4.0         0 DP-1
## 2     1    12.3      13.1 5.11   1 3.9         0 DP-2
## 3     1    12.3      13.1 5.11   1 3.5         0 DP-3
## 4     1    12.3      13.1 5.11   1 3.6         0 DP-4
## 5     1    12.3      13.1 5.11   F 3.9         0 DP-5
## 6     1    12.3      13.1 5.11   F 4.0         0 DP-6
              meta1$swab_time<-with(meta1,paste(swab,visit,sep="_"))
              
              
            #Copy Across Metadata  
              sample2<-sample1 %>% dplyr::select(swab_time) %>% left_join(meta1,"swab_time")
              #head(sample2)
              sample2$site<- sapply(strsplit(sample2$swab,"-"),"[",1)
              rownames(sample2)<-sample1[,1]
              
            ##Sample Type
              sample2$sampletype<-"swab"
              sample2$sampletype[grep("MC",rownames(sample2))]<-"mock"
              sample2$sampletype[grep("NEG",rownames(sample2))]<-"negative_control"
              
            #Sequence Abundance
              otu1<-readRDS('Newt 2021 SeqTab.rds')
              #head(otu1)
              
              #Infection Data
                infection<-read.csv('newt infection codes.csv',header=T)
                head(infection)
##   swab lesion timecode
## 1 DP-1   <NA>        1
## 2 DP-2   <NA>        1
## 3 DP-3   <NA>        1
## 4 DP-4   <NA>        1
## 5 DP-5   <NA>        1
## 6 DP-6   <NA>        1
                infection$swab_time<-with(infection,paste(swab,timecode,sep="_"))
                sample2$infectioncode<-infection$lesion[match(sample2$swab_time,infection$swab_time)]
                  head(sample2)
##       swab_time visit airtemp watertemp   ph sex svl infection  swab site
## DP-1     DP-1_1     1    12.3      13.1 5.11   1 4.0         0  DP-1   DP
## DP-10   DP-10_1     1    12.3      13.1 5.11   1 3.5         0 DP-10   DP
## DP-11   DP-11_1     1    12.3      13.1 5.11   1 3.5         0 DP-11   DP
## DP-12   DP-12_1     1    12.3      13.1 5.11   1 3.3         0 DP-12   DP
## DP-13   DP-13_1     1    12.3      13.1 5.11   1 3.6         0 DP-13   DP
## DP-14   DP-14_1     1    12.3      13.1 5.11   1 3.9         0 DP-14   DP
##       sampletype infectioncode
## DP-1        swab          <NA>
## DP-10       swab          <NA>
## DP-11       swab          <NA>
## DP-12       swab          <NA>
## DP-13       swab          <NA>
## DP-14       swab          <NA>
            #Re-Name Rows with COrrect Sample IDs
              rownames(otu1)<-sample1[,1]

Phyloseq Build

newt <- phyloseq(otu_table(otu1,taxa_are_rows = FALSE),   # create new phyloseq otu table, call it "newt"
                 sample_data(sample2), 
                 tax_table((tax1)))
newt 
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 13903 taxa and 152 samples ]
## sample_data() Sample Data:       [ 152 samples by 12 sample variables ]
## tax_table()   Taxonomy Table:    [ 13903 taxa by 7 taxonomic ranks ]

Decontamination

         ##### REMOVE  MOCK SAMPLES
              
              newt.nomock<-prune_samples(sample_data(newt)$sampletype %in% c("swab","negative_control"),newt)
              newt.nomock
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 13903 taxa and 150 samples ]
## sample_data() Sample Data:       [ 150 samples by 12 sample variables ]
## tax_table()   Taxonomy Table:    [ 13903 taxa by 7 taxonomic ranks ]
              #colnames(tax_table(newt.nomock))<-c("Kingdom","Phylum","Class","Order","Family","Genus","Species")
              
              
        ###### REMOVE NEGATIVE CONTAMINATION
              
              sample_data(newt.nomock)$is.neg <- ifelse(grepl("negative_control",sample_data(newt.nomock)$sampletype),TRUE,FALSE)
              contamdf.prev <- isContaminant(newt.nomock, method="prevalence", neg="is.neg",threshold=0.5)
              table(contamdf.prev$contaminant)
## 
## FALSE  TRUE 
## 13870    33
                # 33 contaminants identified 
              
              ###### Plot of Contaminant Frequency 
                    # Make phyloseq object of presence-absence in negative controls and true samples
                        ps.pa <- transform_sample_counts(newt.nomock, function(abund) 1*(abund>0))
                        ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
                        ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)
                        
                    # Make data.frame of prevalence in positive and negative samples
                        df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                                        contaminant=contamdf.prev$contaminant)
                        
                          ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() +
                            xlab("Prevalence (Negative Controls)") + ylab("Prevalence (True Samples)")

              ###### Clean Contaminants Out
                          newt.clean<-prune_taxa(contamdf.prev$contaminant==FALSE,newt.nomock)
                          ntaxa(newt.nomock) - ntaxa(newt.clean)
## [1] 33

Further Cleaning

        ##### REMOVE NEGATIVES SAMPLES
              
              newt.clean<-prune_samples(sample_data(newt.clean)$sampletype=="swab",newt.clean)
              newt.clean
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 13870 taxa and 148 samples ]
## sample_data() Sample Data:       [ 148 samples by 13 sample variables ]
## tax_table()   Taxonomy Table:    [ 13870 taxa by 7 taxonomic ranks ]
        #Prune Taxa With No Phylum Assignment 
          ps_prune<-prune_taxa(as.vector(!is.na(tax_table(newt.clean)[,2])),newt.clean)
          ntaxa(newt.clean)-ntaxa(ps_prune) #232 lost
## [1] 232
        #Prune Chloroplasts
          ps_prune_nochloro<-prune_taxa(as.vector(tax_table(ps_prune)[,3]!="Chloroplast"),ps_prune)
          ntaxa(ps_prune)-ntaxa(ps_prune_nochloro) #431 lost
## [1] 431
        #Remove Mitochondria    
          ps_prune_nochloro_nomito<-prune_taxa(as.vector(tax_table(ps_prune_nochloro)[,5]!="Mitochondria"),ps_prune_nochloro)
        
        #Remove Archaea 
          ps_prune_nochloro_noarchaea<-prune_taxa(as.vector(tax_table(ps_prune_nochloro_nomito)[,1]!="Archaea"),ps_prune_nochloro_nomito)
            ntaxa(ps_prune_nochloro_noarchaea)-ntaxa(ps_prune_nochloro_nomito) #93 lost
## [1] -93

Filter Based on Prevalence and Abundance

####### Read Threshold and Sample Prevalence Threshold
    ps_clean_filter = filter_taxa(ps_prune_nochloro_noarchaea, function(x) sum(x > 20) > (0.02*length(x)), TRUE)

Mock Community

##Inspect Mock After Removing ASVs filtered at previous step 
  ps_mock<-prune_samples(sample_data(newt)$sampletype %in% c("mock"),newt)
  ps_mock<-filter_taxa(ps_mock, function(x) (sum(x) > 0),TRUE)
  data.frame(reads=taxa_sums(ps_mock),genus=unname(tax_table(ps_mock))[,6])
##                                                                                                                                                                                                                                                               reads
## TACGGAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGCTTGATAAGCCGGTTGTGAAAGCCCCGGGCTCAACCTGGGAACGGCATCCGGAACTGTCAAGCTAGAGTGCAGGAGAGGAAGGTAGAATTCCCGGTGTAGCGGTGAAATGCGTAGAGATCGGGAGGAATACCAGTGGCGAAGGCGGCCTTCTGGACTGACACTGACACTGAGGTGCGAAAGCGTGGGTAGCAAACAGG    11
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG  1322
## TACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTAAGAGGTGAAAGCCCAGAGCTCAACTCTGGAATTGCCTTTTAGACTGCATCGCTTGAATCATGGAGAGGTCAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTGACTGGACATGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGG 21387
## TACAGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGTTCGTTAAGTTGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCAAAACTGTCGAGCTAGAGTATGGTAGAGGGTGGTGGAATTTCCTGTGTAGCGGTGAAATGCGTAGATATAGGAAGGAACACCAGTGGCGAAGGCGACCACCTGGACTGATACTGACACTGAGGTGCGAAAGCGTGGGGAGCAAACAGG 13464
## TACGGAGGATCCAAGCGTTATCCGGAATCATTGGGTTTAAAGGGTCCGTAGGCGGTTTAATAAGTCAGTGGTGAAAGCCCATCGCTCAACGGTGGAACGGCCATTGATACTGTTAGACTTGAATTATTAGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGAGATTACATGGAATACCAATTGCGAAGGCAGGTTACTACTAATTGATTGACGCTGATGGACGAAAGCGTGGGTAGCGAACAGG 15499
## TACGAAGGGGGCTAGCGTTGTTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGATATTTAAGTCAGGGGTGAAATCCCAGAGCTCAACTCTGGAACTGCCTTTGATACTGGGTATCTTGAGTATGGAAGAGGTAAGTGGAATTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAGGAACACCAGTGGCGAAGGCGGCTTACTGGTCCATAACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGG  6853
## TACGGAGGATCCAAGCGTTATCCGGAATCATTGGGTTTAAAGGGTCCGTAGGCGGTTTAATAAGTCAGTGGTGAAAGCCCATCGCTCAACGGTGGAACGGCCATTGATACTGTTAGACTTGAATTATTAGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGAGATTACATGGAATACCAATTGCGAAGGCAGGTTACTACTAGTTGATTGACGCTGATGGACGAAAGCGTGGGTAGCGAACAGG  2507
## TACGTAGGGTCCGAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAATCTGGAGGCTCAACCTCCAGCCTGCAGTGGGTACGGGCAGACTAGAGTGCGGTAGGGGAGAATGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGTTCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACAGG  2248
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCCGAAACTGGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG  1676
## TACGTAGGTGGCAAGCGTTATCCGGAATTATTGGGCGTAAAGCGCGCGCAGGTGGTTTCTTAAGTCTGATGTGAAAGCCCACGGCTCAACCGTGGAGGGTCATTGGAAACTGGGAGACTTGAGTGCAGAAGAGGAAAGTGGAATTCCATGTGTAGCGGTGAAATGCGTAGAGATATGGAGGAACACCAGTGGCGAAGGCGACTTTCTGGTCTGTAACTGACACTGAGGCGCGAAAGCGTGGGGAGCAAACAGG  1790
## TACGGAGGATCCGAGCGTTATCCGGAATTATTGGGTTTAAAGGGTTCGTAGGCGGTCAGATAAGTCAGTGGTGAAATTTCCTAGCTTAACTAGGACACGGCCATTGATACTGTTTGACTTGAATAGTATGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGATATTACATGGAATACCAATTGCGAAGGCAGGTTACTACGTACTTATTGACGCTGATGAACGAAAGCGTGGGGAGCGAACAGG  1364
## TACGGAGGATCCGAGCGTTATCCGGAATTATTGGGTTTAAAGGGTTCGTAGGCGGTCAGATAAGTCAGTGGTGAAATTTCTTAGCTTAACTAAGACACGGCCATTGATACTGTTTGACTTGAATAGTATGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGATATTACATGGAATACCAATTGCGAAGGCAGGTTACTACGTACTTATTGACGCTGATGAACGAAAGCGTGGGGAGCGAACAGG  1085
## TACGTAGGTGGCAAGCGTTATCCGGAATTATTGGGCGTAAAGCGCGCGTAGGCGGTTTCTTAAGTCTGATGTGAAAGCCCACGGCTCAACCGTGGAGGGTCATTGGAAACTGGGAAACTTGAGTGCAGAAGAGGAAAGTGGAATTCCATGTGTAGCGGTGAAATGCGCAGAGATATGGAGGAACACCAGTGGCGAAGGCGACTTTCTGGTCTGTAACTGACGCTGATGTGCGAAAGCGTGGGGATCAAACAGG   893
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGGCAGGCTTGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG   932
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCATGCAGGTGGTTTGTTAAGTCAGATGTGAAAGCCCTGGGCTCAACCTAGGAATCGCATTTGAAACTGACAAGCTAGAGTACTGTAGAGGGGGGTAGAATTTCAGGTGTAGCGGTGAAATGCGTAGAGATCTGAAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAGATACTGACACTCAGATGCGAAAGCGTGGGGAGCAAACAGG   857
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTTGTTAAGTCAGATGTGAAATCCCCGAGCTCAACTTGGGAACTGCATTTGAAACTGGCAAGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG   787
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTTGTTAAGTCAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCTGATACTGGCAAGCTTGAGTCTCGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACGAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG   703
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTGATTAAGTCAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCTGATACTGGTCAGCTTGAGTCTCGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACGAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG   654
## TACGGAGGGTGCGAGCGTTAATCGGAATAACTGGGCGTAAAGGGCACGCAGGCGGTTGATTAAGTGAGATGTGAAAGCCCCGGGCTTAACCTGGGAATTGCATTTCATACTGGTCAACTAGAGTACTTTAGGGAGGGGTAGAATTCCACGTGTAGCGGTGAAATGCGTAGAGATGTGGAGGAATACCGAAGGCGAAGGCAGCCCCTTGGGAATGTACTGACGCTCATGTGCGAAAGCGTGGGGAGCAAACAGG   438
## TACAGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGTTCGTTAAGTTGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCCAAAACTGGCGAGCTAGAGTAGGGCAGAGGGTGGTGGAATTTCCTGTGTAGCGGTGAAATGCGTAGATATAGGAAGGAACACCAGTGGCGAAGGCGACCACCTGGGCTCATACTGACACTGAGGTGCGAAAGCGTGGGGAGCAAACAGG   323
## TACGTAGGGCGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAAGCCCGGGGCTTAACCCCGGGTGTGCAGTGGGTACGGGCAGACTAGAGTGCAGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGGTCTCTGGGCTGTTACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACAGG   200
## TACGGAGGATCCGAGCGTTATCCGGAATTATTGGGTTTAAAGGGTTCGTAGGCGGCTTTGTAAGTCAGTGGTGAAATTTCCTAGCTTAACTAGGACACTGCCATTGATACTGCAGAGCTTGAATAATATGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGATATTACATGGAATACCAATTGCGAAGGCAGGTTACTACGTATTTATTGACGCTGATGAACGAAAGCGTGGGGAGCGAACAGG   278
## TACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTTGATAAGCCAGATGTGAAATCCCCGAGCTCAACTTGGGAACTGCGTTTGGAACTGTCAGACTAGAGTGCGTCAGAGGGGGGTGGAATTCCGCGTGTAGCAGTGAAATGCGTAGAGATGCGGAGGAACACCGATGGCGAAGGCAGCCCCCTGGGATGACACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACAGG    76
## TACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAATCCCGAGGCTCAACCTCGGGTCTGCAGTGGGTACGGGCAGACTAGAGTGCGGTAGGGGAGATTGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGATCTCTGGGCCGCTACTGACGCTGAGGAGCGAAAGGGTGGGGAGCAAACAGG    54
## CGATGAACATGATTAGCGGTACGACTTTGGCCACCGTGGTCAGTTGATTGATCAGCGCCGCCTCCTTGATCCCGCGCATCACCAAAAAATGCACGGCCCACAGCAGCACCGACGCGCAGCCGATGGCGATCGGCGTATTGCCTTGGCCAAACACCGGAAAGAAGTAGCCGAGGGTGCTGAACAGCAGGACGAAATAACCGACGTTGCCCAGCCACGCACTGATCCAGTAACCCCAGGCTG                  6
##                                                                                                                                                                                                                                                                                ta1
## TACGGAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGCTTGATAAGCCGGTTGTGAAAGCCCCGGGCTCAACCTGGGAACGGCATCCGGAACTGTCAAGCTAGAGTGCAGGAGAGGAAGGTAGAATTCCCGGTGTAGCGGTGAAATGCGTAGAGATCGGGAGGAATACCAGTGGCGAAGGCGGCCTTCTGGACTGACACTGACACTGAGGTGCGAAAGCGTGGGTAGCAAACAGG            Halomonas
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG           Klebsiella
## TACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTAAGAGGTGAAAGCCCAGAGCTCAACTCTGGAATTGCCTTTTAGACTGCATCGCTTGAATCATGGAGAGGTCAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTGACTGGACATGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGG         Sphingomonas
## TACAGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGTTCGTTAAGTTGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCAAAACTGTCGAGCTAGAGTATGGTAGAGGGTGGTGGAATTTCCTGTGTAGCGGTGAAATGCGTAGATATAGGAAGGAACACCAGTGGCGAAGGCGACCACCTGGACTGATACTGACACTGAGGTGCGAAAGCGTGGGGAGCAAACAGG          Pseudomonas
## TACGGAGGATCCAAGCGTTATCCGGAATCATTGGGTTTAAAGGGTCCGTAGGCGGTTTAATAAGTCAGTGGTGAAAGCCCATCGCTCAACGGTGGAACGGCCATTGATACTGTTAGACTTGAATTATTAGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGAGATTACATGGAATACCAATTGCGAAGGCAGGTTACTACTAATTGATTGACGCTGATGGACGAAAGCGTGGGTAGCGAACAGG       Flavobacterium
## TACGAAGGGGGCTAGCGTTGTTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGATATTTAAGTCAGGGGTGAAATCCCAGAGCTCAACTCTGGAACTGCCTTTGATACTGGGTATCTTGAGTATGGAAGAGGTAAGTGGAATTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAGGAACACCAGTGGCGAAGGCGGCTTACTGGTCCATAACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGG         Neorhizobium
## TACGGAGGATCCAAGCGTTATCCGGAATCATTGGGTTTAAAGGGTCCGTAGGCGGTTTAATAAGTCAGTGGTGAAAGCCCATCGCTCAACGGTGGAACGGCCATTGATACTGTTAGACTTGAATTATTAGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGAGATTACATGGAATACCAATTGCGAAGGCAGGTTACTACTAGTTGATTGACGCTGATGGACGAAAGCGTGGGTAGCGAACAGG       Flavobacterium
## TACGTAGGGTCCGAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAATCTGGAGGCTCAACCTCCAGCCTGCAGTGGGTACGGGCAGACTAGAGTGCGGTAGGGGAGAATGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGTTCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACAGG            Agromyces
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCCGAAACTGGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG          Citrobacter
## TACGTAGGTGGCAAGCGTTATCCGGAATTATTGGGCGTAAAGCGCGCGCAGGTGGTTTCTTAAGTCTGATGTGAAAGCCCACGGCTCAACCGTGGAGGGTCATTGGAAACTGGGAGACTTGAGTGCAGAAGAGGAAAGTGGAATTCCATGTGTAGCGGTGAAATGCGTAGAGATATGGAGGAACACCAGTGGCGAAGGCGACTTTCTGGTCTGTAACTGACACTGAGGCGCGAAAGCGTGGGGAGCAAACAGG             Bacillus
## TACGGAGGATCCGAGCGTTATCCGGAATTATTGGGTTTAAAGGGTTCGTAGGCGGTCAGATAAGTCAGTGGTGAAATTTCCTAGCTTAACTAGGACACGGCCATTGATACTGTTTGACTTGAATAGTATGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGATATTACATGGAATACCAATTGCGAAGGCAGGTTACTACGTACTTATTGACGCTGATGAACGAAAGCGTGGGGAGCGAACAGG             Myroides
## TACGGAGGATCCGAGCGTTATCCGGAATTATTGGGTTTAAAGGGTTCGTAGGCGGTCAGATAAGTCAGTGGTGAAATTTCTTAGCTTAACTAAGACACGGCCATTGATACTGTTTGACTTGAATAGTATGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGATATTACATGGAATACCAATTGCGAAGGCAGGTTACTACGTACTTATTGACGCTGATGAACGAAAGCGTGGGGAGCGAACAGG             Myroides
## TACGTAGGTGGCAAGCGTTATCCGGAATTATTGGGCGTAAAGCGCGCGTAGGCGGTTTCTTAAGTCTGATGTGAAAGCCCACGGCTCAACCGTGGAGGGTCATTGGAAACTGGGAAACTTGAGTGCAGAAGAGGAAAGTGGAATTCCATGTGTAGCGGTGAAATGCGCAGAGATATGGAGGAACACCAGTGGCGAAGGCGACTTTCTGGTCTGTAACTGACGCTGATGTGCGAAAGCGTGGGGATCAAACAGG       Staphylococcus
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGGCAGGCTTGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG           Salmonella
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCATGCAGGTGGTTTGTTAAGTCAGATGTGAAAGCCCTGGGCTCAACCTAGGAATCGCATTTGAAACTGACAAGCTAGAGTACTGTAGAGGGGGGTAGAATTTCAGGTGTAGCGGTGAAATGCGTAGAGATCTGAAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAGATACTGACACTCAGATGCGAAAGCGTGGGGAGCAAACAGG               Vibrio
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTTGTTAAGTCAGATGTGAAATCCCCGAGCTCAACTTGGGAACTGCATTTGAAACTGGCAAGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG          Plesiomonas
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTTGTTAAGTCAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCTGATACTGGCAAGCTTGAGTCTCGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACGAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG Escherichia-Shigella
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTGATTAAGTCAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCTGATACTGGTCAGCTTGAGTCTCGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACGAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG                 <NA>
## TACGGAGGGTGCGAGCGTTAATCGGAATAACTGGGCGTAAAGGGCACGCAGGCGGTTGATTAAGTGAGATGTGAAAGCCCCGGGCTTAACCTGGGAATTGCATTTCATACTGGTCAACTAGAGTACTTTAGGGAGGGGTAGAATTCCACGTGTAGCGGTGAAATGCGTAGAGATGTGGAGGAATACCGAAGGCGAAGGCAGCCCCTTGGGAATGTACTGACGCTCATGTGCGAAAGCGTGGGGAGCAAACAGG       Actinobacillus
## TACAGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGTTCGTTAAGTTGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCCAAAACTGGCGAGCTAGAGTAGGGCAGAGGGTGGTGGAATTTCCTGTGTAGCGGTGAAATGCGTAGATATAGGAAGGAACACCAGTGGCGAAGGCGACCACCTGGGCTCATACTGACACTGAGGTGCGAAAGCGTGGGGAGCAAACAGG          Pseudomonas
## TACGTAGGGCGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAAGCCCGGGGCTTAACCCCGGGTGTGCAGTGGGTACGGGCAGACTAGAGTGCAGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGGTCTCTGGGCTGTTACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACAGG              Kocuria
## TACGGAGGATCCGAGCGTTATCCGGAATTATTGGGTTTAAAGGGTTCGTAGGCGGCTTTGTAAGTCAGTGGTGAAATTTCCTAGCTTAACTAGGACACTGCCATTGATACTGCAGAGCTTGAATAATATGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGATATTACATGGAATACCAATTGCGAAGGCAGGTTACTACGTATTTATTGACGCTGATGAACGAAAGCGTGGGGAGCGAACAGG             Myroides
## TACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTTGATAAGCCAGATGTGAAATCCCCGAGCTCAACTTGGGAACTGCGTTTGGAACTGTCAGACTAGAGTGCGTCAGAGGGGGGTGGAATTCCGCGTGTAGCAGTGAAATGCGTAGAGATGCGGAGGAACACCGATGGCGAAGGCAGCCCCCTGGGATGACACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACAGG            Vogesella
## TACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAATCCCGAGGCTCAACCTCGGGTCTGCAGTGGGTACGGGCAGACTAGAGTGCGGTAGGGGAGATTGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGATCTCTGGGCCGCTACTGACGCTGAGGAGCGAAAGGGTGGGGAGCAAACAGG         Plantibacter
## CGATGAACATGATTAGCGGTACGACTTTGGCCACCGTGGTCAGTTGATTGATCAGCGCCGCCTCCTTGATCCCGCGCATCACCAAAAAATGCACGGCCCACAGCAGCACCGACGCGCAGCCGATGGCGATCGGCGTATTGCCTTGGCCAAACACCGGAAAGAAGTAGCCGAGGGTGCTGAACAGCAGGACGAAATAACCGACGTTGCCCAGCCACGCACTGATCCAGTAACCCCAGGCTG                              <NA>

##Library Sizes After Cleaning

#Summary
 range(sample_sums(ps_clean_filter)) 
## [1]  1278 69644
  mean(sample_sums(ps_clean_filter)) 
## [1] 22264.95
#Table
  sample_sums(ps_clean_filter)
##    DP-1   DP-10   DP-11   DP-12   DP-13   DP-14   DP-15   DP-16   DP-17   DP-18 
##   28480   38000   33723   36636   56889   42389   32158   34777   27873   28435 
##   DP-19   DP-1C    DP-2   DP-20   DP-21   DP-22   DP-23   DP-24   DP-25   DP-26 
##   32201   14498   43533   32864   39365   36527   31885   42616   52852   43414 
##   DP-27   DP-2C    DP-3   DP-3C    DP-4   DP-4C    DP-5   DP-5C    DP-6    DP-7 
##   27425   25949   35236   17054   69544   20893   53744   15085   38979   55027 
##    DP-8   DP-8C    DP-9   R10-1  R10-10  R10-11  R10-12  R10-13  R10-14  R10-15 
##   69644   14214   24015   11275    8525   16545   14457    9632   17585   14677 
##  R10-16  R10-17  R10-18  R10-19  R10-1C   R10-2  R10-20  R10-2C   R10-3  R10-3C 
##   19851   12916   11771   11105   15820   15172   15998   28945    9717   28280 
##   R10-4   R10-5   R10-6  R10-6C   R10-7  R10-7C   R10-8  R10-8C   R4-10   R4-11 
##    9688   10696   13754    9088    9817   18480   13746   19026    7882   26514 
##   R4-12  R4-12C   R4-13  R4-13C   R4-14  R4-14C  R4-18C   R4-1C    R4-2   R4-21 
##   10092   21847   20223   17351   15567   23889    9373   11501   16158   27400 
##  R4-21C   R4-22  R4-22C   R4-23  R4-23C   R4-24   R4-25   R4-26   R4-27   R4-28 
##   21423    7133   16220   19078   18841    9133   12559   13181    6693   19605 
##   R4-2C    R4-3   R4-30   R4-31   R4-32   R4-33   R4-34   R4-35   R4-36   R4-38 
##    8476   28765    9835    7668    1278   20874   13017    7076   29156   24914 
##   R4-39   R4-3C    R4-4   R4-40   R4-41   R4-42   R4-43   R4-44   R4-45  R4-45C 
##   22101   13156   51135   14121   11672   10043   18843   31806   10400   27929 
##   R4-46   R4-47   R4-48   R4-49   R4-4C   R4-50   R4-5C    R4-6   R4-6C    R4-7 
##   22375   45324    5989   24170   18693   24692   30438   27978   16159   21036 
##    R4-8    R4-9   R44-1  Rnew-1 Rnew-10 Rnew-11 Rnew-12 Rnew-13 Rnew-14 Rnew-15 
##   15286   16207   11206   17843   18919   19913   28563   21616   23889   36707 
##  Rnew-2  Rnew-3 Rnew-37 Rnew-38 Rnew-39  Rnew-4 Rnew-40 Rnew-41 Rnew-42 Rnew-43 
##   28399   41446   23936   23781    8069   27542    8690   21795   25391   13301 
## Rnew-44 Rnew-45 Rnew-46 Rnew-47 Rnew-48 Rnew-49  Rnew-5 Rnew-50 Rnew-51 Rnew-52 
##   21870   14670    5943    8300    7312   12217   49985   11141    5390    8642 
## Rnew-53 Rnew-54 Rnew-55 Rnew-56  Rnew-6  Rnew-7  Rnew-8  Rnew-9 
##    7532    8148   20280   19401   38318   18771   52717   38805

Prevalence Models

Probability of Observing C and D Infection

sample2$C_or_D<-ifelse(sample2$infectioncode %in% c("C","D"),1,0)
#table(sample2$infectioncode)
with(sample2,table(infectioncode,site))
##              site
## infectioncode DP R10 R4 Rnew
##             A  0   0 11   10
##             B  0   0 14    6
##             C  0   0  9    4
##             D  0   0  2    0
prev_mod1<-brm(C_or_D ~ site,data = sample2,family=bernoulli()) 
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.12 seconds (Warm-up)
## Chain 1:                0.126 seconds (Sampling)
## Chain 1:                0.246 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 2e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.124 seconds (Warm-up)
## Chain 2:                0.102 seconds (Sampling)
## Chain 2:                0.226 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 4e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.109 seconds (Warm-up)
## Chain 3:                0.146 seconds (Sampling)
## Chain 3:                0.255 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 2e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.094 seconds (Warm-up)
## Chain 4:                0.109 seconds (Sampling)
## Chain 4:                0.203 seconds (Total)
## Chain 4:
## Warning: There were 2 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess
  summary(prev_mod1)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: C_or_D ~ site 
##    Data: sample2 (Number of observations: 147) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Regression Coefficients:
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept   -21.33     34.53   -92.28    -3.67 1.02      145      133
## siteR10     -10.76     67.42  -128.74    43.14 1.02      232      126
## siteR4       19.95     34.53     2.27    91.20 1.02      144      134
## siteRnew     19.20     34.53     1.46    90.16 1.02      146      133
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

ALPHA DIVERSITY

Alpha Calculations First Visit

##Rarefy 
  newt1<-prune_samples(sample_data(ps_clean_filter)$visit==1 & sample_data(ps_clean_filter)$site!="R44",ps_clean_filter)
  newt1 #121 samples after removing R44
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 893 taxa and 121 samples ]
## sample_data() Sample Data:       [ 121 samples by 13 sample variables ]
## tax_table()   Taxonomy Table:    [ 893 taxa by 7 taxonomic ranks ]
#Rarefy data 5000 reads 
  newt1r<-rarefy_even_depth(newt1,rngseed = 01042020,sample.size = 5000) 
## `set.seed(1042020)` was used to initialize repeatable random subsampling.
## Please record this for your records so others can reproduce.
## Try `set.seed(1042020); .Random.seed` for the full vector
## ...
## 1 samples removedbecause they contained fewer reads than `sample.size`.
## Up to first five removed samples are:
## R4-32
## ...
## 6OTUs were removed because they are no longer 
## present in any sample after random subsampling
## ...
  #range(sample_sums(newt1r))

##Estimate alpha diversity metrics
    newt1_rich<-estimate_richness(newt1r,measures=c("Observed"))
    newt1_rich$swab<-gsub("\\.","-",rownames(newt1_rich)) 
    newt1_rich$swab_time<-with(newt1_rich,paste(swab, "1",sep="_"))
##Join  
  newt1_rich_meta<-left_join(newt1_rich,sample2,"swab_time")
  head(newt1_rich_meta)
##   Observed swab.x swab_time visit airtemp watertemp   ph sex svl infection
## 1       81   DP-1    DP-1_1     1    12.3      13.1 5.11   1 4.0         0
## 2      112  DP-10   DP-10_1     1    12.3      13.1 5.11   1 3.5         0
## 3      143  DP-11   DP-11_1     1    12.3      13.1 5.11   1 3.5         0
## 4       55  DP-12   DP-12_1     1    12.3      13.1 5.11   1 3.3         0
## 5      111  DP-13   DP-13_1     1    12.3      13.1 5.11   1 3.6         0
## 6       88  DP-14   DP-14_1     1    12.3      13.1 5.11   1 3.9         0
##   swab.y site sampletype infectioncode C_or_D
## 1   DP-1   DP       swab          <NA>      0
## 2  DP-10   DP       swab          <NA>      0
## 3  DP-11   DP       swab          <NA>      0
## 4  DP-12   DP       swab          <NA>      0
## 5  DP-13   DP       swab          <NA>      0
## 6  DP-14   DP       swab          <NA>      0
##Standardise and Factorise 
  newt1_rich_meta$z_svl<- with(newt1_rich_meta,as.numeric(scale(svl))) # mean centered svl (snout-ventral length)
  newt1_rich_meta$watertemp<-factor(newt1_rich_meta$watertemp)
  newt1_rich_meta$infection<-factor(newt1_rich_meta$infection)
  #newt1_rich_meta$ph<-factor(newt1_rich_meta$ph)
  newt1_rich_meta$site<-factor(newt1_rich_meta$site)

Alpha Models

Site, Sex, SVL

#Maximal Model 
  rich_mod1<-brm(Observed ~  site + sex + z_svl, data = newt1_rich_meta, family = negbinomial())
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 3.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.07 seconds (Warm-up)
## Chain 1:                0.074 seconds (Sampling)
## Chain 1:                0.144 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 7e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.071 seconds (Warm-up)
## Chain 2:                0.068 seconds (Sampling)
## Chain 2:                0.139 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 7e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.071 seconds (Warm-up)
## Chain 3:                0.068 seconds (Sampling)
## Chain 3:                0.139 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 6e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.068 seconds (Warm-up)
## Chain 4:                0.067 seconds (Sampling)
## Chain 4:                0.135 seconds (Total)
## Chain 4:
  summary(rich_mod1)
##  Family: negbinomial 
##   Links: mu = log 
## Formula: Observed ~ site + sex + z_svl 
##    Data: newt1_rich_meta (Number of observations: 120) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Regression Coefficients:
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     4.58      0.07     4.44     4.72 1.00     3517     2815
## siteR10       0.70      0.14     0.44     0.97 1.00     2573     2392
## siteR4       -0.55      0.09    -0.72    -0.38 1.00     3011     2881
## siteRnew     -0.74      0.09    -0.92    -0.55 1.00     2801     2769
## sexF         -0.02      0.09    -0.20     0.16 1.00     2992     2249
## sexM         -0.19      0.18    -0.55     0.15 1.00     2673     2437
## z_svl        -0.02      0.04    -0.10     0.07 1.00     3060     2836
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape     9.53      1.44     6.96    12.57 1.00     3495     3207
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#Diagnostics
  pp_check(rich_mod1)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

  bayes_R2(rich_mod1)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.7427242 0.03864956 0.6462914 0.7939698
#Plot
  conditional_effects(rich_mod1)
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"

## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"

## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"

#Null Model 
    rich_mod_null<-brm(Observed ~  1, data = newt1_rich_meta, family = negbinomial())
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 2.2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.045 seconds (Warm-up)
## Chain 1:                0.041 seconds (Sampling)
## Chain 1:                0.086 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 9e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.044 seconds (Warm-up)
## Chain 2:                0.047 seconds (Sampling)
## Chain 2:                0.091 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 9e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.045 seconds (Warm-up)
## Chain 3:                0.051 seconds (Sampling)
## Chain 3:                0.096 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 9e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.044 seconds (Warm-up)
## Chain 4:                0.047 seconds (Sampling)
## Chain 4:                0.091 seconds (Total)
## Chain 4:
#Add ICs
    rich_mod1<-add_criterion(rich_mod1,"loo")
    rich_mod_null<-add_criterion(rich_mod_null,"loo")
#Compare    
    loo_compare(rich_mod1,rich_mod_null)
##               elpd_diff se_diff
## rich_mod1       0.0       0.0  
## rich_mod_null -65.3       9.4
#Extract Estimates
    rich_predict<-  conditional_effects(rich_mod1)[[1]]

pH

  #Numeric
    #newt1_rich_meta$ph<-as.numeric(as.character(newt1_rich_meta$ph))

  #Model 
    ph_model1<-brm(Observed ~ ph + I(ph^2) + infection,data=newt1_rich_meta,family=negbinomial())
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 3.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.739 seconds (Warm-up)
## Chain 1:                0.695 seconds (Sampling)
## Chain 1:                1.434 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 7e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.659 seconds (Warm-up)
## Chain 2:                0.642 seconds (Sampling)
## Chain 2:                1.301 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 6e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.669 seconds (Warm-up)
## Chain 3:                0.65 seconds (Sampling)
## Chain 3:                1.319 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 6e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.689 seconds (Warm-up)
## Chain 4:                0.791 seconds (Sampling)
## Chain 4:                1.48 seconds (Total)
## Chain 4:
    summary(ph_model1)
##  Family: negbinomial 
##   Links: mu = log 
## Formula: Observed ~ ph + I(ph^2) + infection 
##    Data: newt1_rich_meta (Number of observations: 120) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Regression Coefficients:
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     41.14      4.08    33.15    49.28 1.00     1724     1834
## ph           -12.89      1.47   -15.82    -9.98 1.00     1732     1883
## IphE2          1.12      0.13     0.86     1.38 1.00     1746     1841
## infection1    -0.24      0.09    -0.41    -0.08 1.00     2491     2341
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape     8.72      1.31     6.36    11.46 1.00     3065     1986
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
    #Diagnostics
      pp_check(ph_model1)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

      bayes_R2(ph_model1)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.7385453 0.03498745 0.6452589 0.7784395
    #Plots
      conditional_effects(ph_model1)
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"

## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"

    #Model Selection 
      ph_model_null<-brm(Observed ~ 1,data=newt1_rich_meta,family=negbinomial())
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 2.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.21 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.043 seconds (Warm-up)
## Chain 1:                0.046 seconds (Sampling)
## Chain 1:                0.089 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 9e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.044 seconds (Warm-up)
## Chain 2:                0.038 seconds (Sampling)
## Chain 2:                0.082 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 8e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.045 seconds (Warm-up)
## Chain 3:                0.043 seconds (Sampling)
## Chain 3:                0.088 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 8e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.047 seconds (Warm-up)
## Chain 4:                0.048 seconds (Sampling)
## Chain 4:                0.095 seconds (Total)
## Chain 4:
      ph_model1<-add_criterion(ph_model1,"loo")
      ph_model_null<-add_criterion(ph_model_null,"loo")
      loo_compare(ph_model_null,ph_model1)
##               elpd_diff se_diff
## ph_model1       0.0       0.0  
## ph_model_null -62.2       8.7
    #Extract   
      ph_predict<-conditional_effects(ph_model1)[[1]]

Infection Class

#Subset Data
  newt1_rich_infection<- newt1_rich_meta %>% mutate(infection=factor(infection)) %>% filter(site %in% c("R4","Rnew"))

#Model 
  infect_mod1<-brm(Observed ~ infection*site,family=negbinomial(), data=newt1_rich_infection)
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 2.6e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.26 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.044 seconds (Warm-up)
## Chain 1:                0.041 seconds (Sampling)
## Chain 1:                0.085 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 6e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.044 seconds (Warm-up)
## Chain 2:                0.034 seconds (Sampling)
## Chain 2:                0.078 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 4e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 2.788 seconds (Warm-up)
## Chain 3:                3.151 seconds (Sampling)
## Chain 3:                5.939 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 4e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.044 seconds (Warm-up)
## Chain 4:                0.037 seconds (Sampling)
## Chain 4:                0.081 seconds (Total)
## Chain 4:
## Warning: There were 133 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: There were 867 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: The largest R-hat is 1.59, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess
    summary(infect_mod1)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Warning: There were 133 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: negbinomial 
##   Links: mu = log 
## Formula: Observed ~ infection * site 
##    Data: newt1_rich_infection (Number of observations: 74) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Regression Coefficients:
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept             -96.94    175.23  -400.40     4.42 1.53        7       25
## infection1            -10.78     18.16   -42.23    -0.06 1.56        7       14
## siteRnew               27.06     47.63    -0.71   109.55 1.55        7       11
## infection1:siteRnew    20.19     34.34     0.02    79.66 1.59        7       13
## 
## Further Distributional Parameters:
##                                                                                     Estimate
## shape 82860327994598638985802098729997548010686563506201240169762507419768653154650947584.00
##                                                                                     Est.Error
## shape 143536241170127931064739043427420160456027114437337950929267711863194371876355833856.00
##       l-95% CI
## shape     5.58
##                                                                                      u-95% CI
## shape 331441311978395688260968415246861568057379908849293253439115776383118818970234257408.00
##       Rhat Bulk_ESS Tail_ESS
## shape 1.52        7       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
    pp_check(infect_mod1)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

    bayes_R2(infect_mod1)
##     Estimate Est.Error          Q2.5     Q97.5
## R2 0.1111109 0.0894021 3.163101e-224 0.2911242
    conditional_effects(infect_mod1)
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"

## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"

  ## Null Model 
          infect_mod_null<-brm(Observed ~ 1,family=negbinomial(), data=newt1_rich_infection)
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 1.8e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.033 seconds (Warm-up)
## Chain 1:                0.031 seconds (Sampling)
## Chain 1:                0.064 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 6e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.032 seconds (Warm-up)
## Chain 2:                0.035 seconds (Sampling)
## Chain 2:                0.067 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 6e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.048 seconds (Warm-up)
## Chain 3:                0.031 seconds (Sampling)
## Chain 3:                0.079 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 6e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.033 seconds (Warm-up)
## Chain 4:                0.027 seconds (Sampling)
## Chain 4:                0.06 seconds (Total)
## Chain 4:
  #Add ICs
      infect_mod1<-add_criterion(infect_mod1,"loo")  
        infect_mod_null<-add_criterion(infect_mod_null,"loo")   
#Select
        print(loo_compare(infect_mod1,infect_mod_null),simplify=F)
##                 elpd_diff  se_diff    elpd_loo   se_elpd_loo p_loo     
## infect_mod_null        0.0        0.0     -325.5        7.8         2.1
## infect_mod1     -1370526.2    85788.1 -1370851.8    85794.0   1370512.0
##                 se_p_loo   looic      se_looic  
## infect_mod_null        0.6      651.0       15.6
## infect_mod1        85789.7  2741703.5   171588.0
#Model Predictions 
    infect_pred_data<-data.frame(site=rep(c("R4","Rnew"),2),infection=rep(c(0,1),each=2))
    infect_predict<-conditional_effects(infect_mod1,"site:infection")[[1]]

Infection Severity

#combine infection categories c/d
  newt1_rich_meta$infection_collapse<-as.factor(newt1_rich_meta$infectioncode)
  table(newt1_rich_meta$infection_collapse)
## 
##  A  B  C  D 
## 19 15 11  2
  levels(newt1_rich_meta$infection_collapse)<-c("A","B","C","C")

#Filter
  newt1_rich_meta_filter<- filter(newt1_rich_meta,!is.na(infection_collapse))
  nrow(newt1_rich_meta_filter) #47
## [1] 47
#Infection Model 
  severity_model1<-brm(Observed ~ infection_collapse,data= newt1_rich_meta_filter,family=negbinomial())
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 2.5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.25 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.024 seconds (Warm-up)
## Chain 1:                0.023 seconds (Sampling)
## Chain 1:                0.047 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 4e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.022 seconds (Warm-up)
## Chain 2:                0.021 seconds (Sampling)
## Chain 2:                0.043 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 4e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.023 seconds (Warm-up)
## Chain 3:                0.022 seconds (Sampling)
## Chain 3:                0.045 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 3e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.022 seconds (Warm-up)
## Chain 4:                0.021 seconds (Sampling)
## Chain 4:                0.043 seconds (Total)
## Chain 4:
  summary(severity_model1)
##  Family: negbinomial 
##   Links: mu = log 
## Formula: Observed ~ infection_collapse 
##    Data: newt1_rich_meta_filter (Number of observations: 47) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Regression Coefficients:
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept               3.96      0.09     3.78     4.14 1.00     3998     2518
## infection_collapseB    -0.03      0.14    -0.29     0.25 1.00     3815     2831
## infection_collapseC    -0.23      0.14    -0.52     0.06 1.00     4226     2941
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape     7.68      1.83     4.72    11.82 1.00     3306     2677
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
    #Diagnostics
    pp_check(severity_model1)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

    bayes_R2(severity_model1)
##      Estimate  Est.Error        Q2.5     Q97.5
## R2 0.07962835 0.05629761 0.003244018 0.2120769
  #Plot
    conditional_effects(severity_model1)
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"

Alpha Plots

Richness by Site

##Richness By Site 
rich_site1 <- ggplot()  + geom_jitter(data=newt1_rich_meta,aes(x=site,y=Observed, fill = site),color="white",shape=21,size=5,width=0.2) + scale_fill_brewer(palette = "Paired") + theme_bw(base_size=15) + theme(text = element_text(size = 18)) #+ guides(fill="none") 
rich_site2<-rich_site1 + geom_errorbar(data=rich_predict,aes(x=site,ymin=lower__,ymax=upper__),width=0.1,lwd=1.2) + geom_point(data=rich_predict,aes(x=site,y=estimate__),shape=23,size=5,fill="white") + labs(x="Site",y="Alpha Diversity \n (Richness)")  + scale_fill_manual(values=sitecols)
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
rich_site3<- rich_site2 +
  theme(legend.key = element_rect(fill = "white"), legend.text = element_text(color = "white"), legend.title = element_text(color = "white")) + guides(fill = guide_legend(override.aes = list(fill = NA)))
rich_site3

ggsave2('Richness by Site Dec23.pdf',rich_site3,height=15,width=16,units="cm")

Shannon Split by Infection Status

rich_infection<- ggplot() + geom_point(data=newt1_rich_infection,aes(x=site,y=Observed,fill=site,shape=infection),position = position_jitterdodge(),size=5,color="white")  + theme_bw(base_size=15) + theme(legend.position = "right") + scale_shape_manual(name="Amphibiothecum \n Infection",values=c(21,24),labels=c("No Visible \n Infection","Infected"))  

rich_infection2<- rich_infection  + labs(y="Alpha Diversity \n (Richness)",x="Site") #+ theme(legend.text = element_text(size=12),legend.title = element_text(size=12))
rich_infection3<- rich_infection2 + geom_errorbar(data=infect_predict,aes(x=site,ymin=lower__,ymax=upper__,group=infection),position = position_dodge(width=0.7),width=0.1,lwd=1.2) + geom_point(data=infect_predict,aes(x=site,y=estimate__,group=infection),position = position_dodge(width=0.7),shape=23,size=5,fill="white") + scale_fill_manual(values=sitecols[c(3,4)]) + guides(fill="none") + guides(shape = guide_legend(override.aes = list(color = "black") ) )+ theme(legend.position = "right")
rich_infection3 

ggsave2('Richness by Infection Dec23.pdf',rich_infection3,height=15,width=18,units="cm")

pH scatterplot

newt1_rich_meta<-mutate(newt1_rich_meta,infection=factor(infection))

richness_ph_scatter<-ggplot(data=newt1_rich_meta, aes(x=as.numeric(as.character(ph)), y=Observed)) +
  geom_jitter(aes(fill=site,shape=infection),color="white", size=5.5, width = 0.05) +
  theme(text = element_text(size = 40)) +
  theme_bw(base_size=15) +
  xlab("pH") + 
  ylab("Alpha Diversity \n (Richness)") + scale_shape_manual(name="Amphibiothecum \n Infection",values=c(21,24),labels=c("No Visible \n Infection","Infected")) + scale_fill_manual(values=sitecols) + labs(fill="Site")

richness_ph_scatter2<- richness_ph_scatter + theme(legend.position = "right")  + scale_x_continuous(n.breaks=8)
richness_ph_scatter3<- richness_ph_scatter2 + geom_ribbon(data=ph_predict,aes(x=ph,ymin=lower__,ymax=upper__),alpha=0.3) + geom_line(data=ph_predict,aes(x=ph,y=estimate__),col="blue") + guides(shape = guide_legend(override.aes = list(color = "black"))) + guides(fill = guide_legend(override.aes = list(shape=21))) + theme(axis.text.x = element_text(angle=45,hjust=1,vjust=1))
richness_ph_scatter3

ggsave2('Richness ph Dec23.pdf',richness_ph_scatter3,height=12,width=20,units="cm")

Combined Alpha Diversity Plot

#Column Format
  alpha_column1<-plot_grid(richness_ph_scatter3,rich_infection3,ncol=1,labels="AUTO",label_size = 30)

ggsave2('Fig_3_Newt Richness Combined Plot.pdf',alpha_column1,width=22,height=40,units="cm")
ggsave2('Fig_3_Newt Richness Combined Plot.tiff',alpha_column1,width=22,height=30,units="cm")

BETA DIVERSITY

##Vegan function to convert phyloseq obj. into something Vegan can call directly
    vegan_otu <- function(physeq) {
      OTU <- otu_table(physeq)
      if (taxa_are_rows(OTU)) {
        OTU <- t(OTU)
      }
      return(as(OTU, "matrix"))
    }

####TIME POINT 1 PCA####

##CLR Transform - Round 1 sampling only
  newt1_clr <- microbiome::transform(newt1, "clr") 
  #head(otu_table(newt1_clr))

##Extract Matrix and Sample Data
  newt1_clr_v<-vegan_otu(newt1_clr) #???
  newt1_clr_s<-as(sample_data(newt1_clr),"data.frame")

##Merge C/D infection codes
    newt1_clr_s$infectioncode[is.na(newt1_clr_s$infectioncode)] <- "0"
    newt1_clr_s$infectioncode<-as.factor(newt1_clr_s$infectioncode)
    levels(newt1_clr_s$infectioncode)<-c("0","A","B","C/D","C/D")
    levels(newt1_clr_s$infectioncode)
## [1] "0"   "A"   "B"   "C/D"
###Principal Components Analysis 
  newt1_pc<-prcomp(newt1_clr_v)
  #biplot(newt1_pc)
  #plot(newt1r_pc)
  summary(newt1_pc)$importance[,1:2]
##                             PC1      PC2
## Standard deviation     15.65805 14.27488
## Proportion of Variance  0.18631  0.15485
## Cumulative Proportion   0.18631  0.34116
  newt1_pc_scores<-scores(newt1_pc)
  newt1_pc_scores_sub<-newt1_pc_scores[,1:2]
  newt1_pc_scores_sub<-cbind(newt1_pc_scores_sub,newt1_clr_s)

##Housekeeping + Label Setup    
  newt1_pc_scores_sub$site<-as.factor(newt1_pc_scores_sub$site)
  newt1_pc_scores_sub$infectioncode<-as.factor(newt1_pc_scores_sub$infectioncode)

### Plot  
  newt1plotdata<-gg_ordiplot(newt1_pc,groups=newt1_clr_s$site,spider=T,plot=T)

  newt1plotdata2<-newt1plotdata$df_ord
  newt1plotdata3<-cbind(newt1plotdata2,newt1_clr_s)
  newt1plotdata3$infection<-factor(newt1plotdata3$infection)
  
  gg1<-ggplot(newt1plotdata3,aes(x=x,y=y)) + geom_segment(data = newt1plotdata$df_spiders, aes(x = cntr.x, xend = x, y = cntr.y, yend = y, color = Group), 
    show.legend = FALSE) + geom_point(aes(shape=infection,fill=site),size=7,color="white") + scale_shape_manual(values=c(21,24),name="Amphibiothecum Infection",labels=c("No Visible Infection","Infected"))
  gg2<- gg1 + geom_path(data = newt1plotdata$df_ellipse, aes(x = x, y = y, color = Group), show.legend = FALSE)
  gg3<- gg2 + theme_bw(base_size=15) + theme(legend.position = "right")
  gg4<- gg3 + guides(fill=guide_legend(nrow=3,byrow=TRUE),shape = guide_legend(override.aes = list(color = "black") ),fill= guide_legend(override.aes = list(color = "black",shape=21))) + labs(x="PC1 (18.6%)",y="PC2 (15.8%)",fill="Site") + guides(fill = guide_legend(override.aes = list(shape=21))) 
## Warning: Duplicated aesthetics after name standardisation: fill
  gg5 <- gg4 + scale_fill_manual(values=sitecols) + scale_color_manual(values=sitecols) 
  gg5

  ggsave2('Newt Visit1 CLR PCA.pdf',gg5,width=20,height=15,units="cm")

PERMANOVA

##Run PERMANOVA on CLR transformed data
raw_permanova<- adonis2(newt1_clr_v ~ site + infection,data=newt1_clr_s,nperm=999,method="euclidean",by="mar")
raw_permanova #site: r2=26.6% , infection: r2=2.8%  
## Permutation test for adonis under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = newt1_clr_v ~ site + infection, data = newt1_clr_s, method = "euclidean", by = "mar", nperm = 999)
##            Df SumOfSqs      R2       F Pr(>F)    
## site        3    42012 0.26605 16.0881  0.001 ***
## infection   1     4449 0.02817  5.1112  0.001 ***
## Residual  116   100974 0.63943                   
## Total     120   157914 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Stacked Barplots

##Filter Samples With No Data using phyloseq object containing clean dataset, which has been filtered from negative/positive control + contaminants, but not filtered for low freq reads.
  # ps_clean_filter_filter<-prune_samples(sample_sums(ps_clean_filter)>0,ps_clean_filter)
  # newt_phylum <- ps_clean_filter_filter %>% aggregate_taxa(level = "Phylum") %>% microbiome::transform(transform = "compositional")
  
#Aggregate To Order Level and Transform to Relative Abundance (fraction of reads, rather than count of reads)
  #newt_order <- newt.clean.bacteria.sub %>% aggregate_taxa(level = "Order") %>% microbiome::transform(transform = "compositional")

#Phylum Colours
  
##Plot Again But With Top 'N' taxa
  newt_phylumcollapse<- newt1 %>% aggregate_taxa(level="Phylum") 
  newt_top5phyla = names(sort(taxa_sums(newt_phylumcollapse), TRUE)[1:5]) #What Are the Names of the most abundant phyla?  
  newt_top5phylum_filter<-subset_taxa(newt1,Phylum %in% newt_top5phyla) #Subset the phyloseq object to those phyla
  infectionlong<-ifelse(sample_data(newt_top5phylum_filter)$infection==0,"Uninfected","Infected") #Site and Infection status- Phylum
  sample_data(newt_top5phylum_filter)$site_infection<-with(sample_data(newt_top5phylum_filter),paste(site,infectionlong,sep=" "))

##Remake our graph but with grouping by site
newt_top5phylum_plot <- newt_top5phylum_filter %>%                   
  aggregate_taxa(level = "Phylum") %>%  
  microbiome::transform(transform = "compositional") %>%
  plot_composition(sample.sort = "Firmicutes", average_by = "site_infection",group_by = "site")
newt_top5phylum_plot <- newt_top5phylum_plot + theme_bw() + 
  scale_fill_manual("Phylum", values=moma.colors("Liu",5)) + 
  theme(legend.position = "bottom") + 
  theme(text = element_text(size = 22)) + labs(x="Site & Infection Status") +
  theme(axis.text.x = element_text(angle=45, vjust=1, hjust=1)) + guides(fill=guide_legend(nrow=3,byrow=TRUE))
newt_top5phylum_plot 

ggsave2('Time 1 Phylum Barplot.pdf',newt_top5phylum_plot,height=20,width=20,units="cm")

CCA Model

Model

    ######
        #Extract Matrix and Sample Data  
          newt1_rare_v<-vegan_otu(newt1)
          newt1_rare_s<-as(sample_data(newt1),"data.frame")

     ##### INDIVIDUAL PREDICTORS
          newt_cca<-cca(newt1_rare_v ~ infection + ph + svl,data=newt1_rare_s,strata=frog_rare_s$frog_id)
            car::Anova(newt_cca)
## Warning in Anova.default(newt_cca): there are rows/columns in vcov.
## that are not in the model matrix:
## CCA1:infection, CCA1:ph, CCA1:svl, CCA2:infection, CCA2:ph, CCA2:svl, CCA3:infection, CCA3:ph, CCA3:svl
## tests may be incorrect
## Analysis of Deviance Table (Type II tests)
## 
## Response: newt1_rare_v
##           Df    Chisq Pr(>Chisq)    
## infection  1  77.0067     <2e-16 ***
## ph         1 446.2242     <2e-16 ***
## svl        1   0.0228       0.88    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

CCA PLot

#Extract Biplot Dat aand Filter Out the Social Group
  bip <- data.frame(scores(newt_cca, choices = 1:2, display = "bp")) 
  bip_filter<- bip  #(bip[-grep("site",rownames(bip)),])

#Scaling Factor
(mul <- ordiArrowMul(bip,display="sites", fill = 0.8))
## [1] 0.8298887
    bip.scl <- bip_filter * 2    
    labs<-c("Infection","pH","SVL")
    
#Extract Plotting Data
  newt_cca_plotdata<- ggvegan::fortify(newt_cca) %>% filter(score=="sites")
  
#Copy Across Metadata
  newt_cca_plotdata$site_id<-newt1_rare_s$site
  newt_cca_plotdata$infection<-ifelse(newt1_rare_s$infection==0,"No Visible Infection","Infected")

#Plot    
  newt_cca_plot1<- ggplot(newt_cca_plotdata,aes(x=CCA1,y=CCA2)) + geom_point(color="white",aes(shape=infection,fill=site_id),size=7) + theme_bw(base_size=15) + scale_shape_manual(values=c(24,21)) + guides(fill = guide_legend(override.aes = list(shape = 21)))   + scale_fill_manual(values=sitecols) + scale_color_manual(values=sitecols) 
#Add Biplot Arrows
  newt_cca_plot2<- newt_cca_plot1 + geom_segment(data=bip.scl,aes(x=0,y=0,xend=CCA1,yend=CCA2),arrow=arrow(length=unit(0.2,"cm")),color="grey") #+ lims(x=c(-6,5))
  newt_cca_plot3 <- newt_cca_plot2 + geom_text_repel(data=bip.scl,aes(x=CCA1,y=CCA2,label=labs),color="grey40",size=5) + labs(fill="Site",shape="Amphibiothecum Infection")  + guides(shape = guide_legend(override.aes = list(color = "black") ))
  newt_cca_plot3

  ggsave2('Newt Visit1 CCA.pdf',newt_cca_plot3,width=22,height=15,units="cm")

Combined Plot

betaplot1<-plot_grid(gg5,newt_cca_plot3,nrow=2,labels="AUTO",rel_widths=c(1,1),label_size = 30)
betaplot1

ggsave2('Newt Visit1 CLR CCA Combined.pdf',betaplot1,width=25,height=25,units="cm")

Core Microbiome Visit 1

Newt Site Specific Core

#Unique Sites
  newtsites <- unique(as.character(meta(newt1)$site))
#Taxonomy Dataframe
  newt1_taxdata<-as.data.frame(tax_table(newt1))

#Write a for loop to go through each of the disease_states one by one and combine identified core taxa into a list.

      sites_core <- c() # an empty object to store information
      sites_core_taxonomy<-c()
      
      for (n in newtsites){ # for each variable n in DiseaseState
          #print(paste0("Identifying Core Taxa for ", n))
          
          ps.sub <- subset_samples(newt1, site == n) # Choose sample from n
          
          core_m <- core_members(ps.sub, # ps.sub is phyloseq selected with only samples from g 
                                 prevalence = 0.9)
          print(paste0("No. of core taxa in ", n, " : ", length(core_m))) # print core taxa identified in each DiseaseState.
          sites_core[[n]] <- core_m # add to a list core taxa for each group.
          sites_core_taxonomy[[n]]<-subset(newt1_taxdata,rownames(newt1_taxdata) %in% core_m)
          #print(list_core)
      }
## [1] "No. of core taxa in DP : 27"
## [1] "No. of core taxa in R10 : 10"
## [1] "No. of core taxa in R4 : 6"
## [1] "No. of core taxa in Rnew : 3"
#Check Core Taxa Not a FUnction of Sample size
      site_ncore_taxa<-data.frame(ntaxa=sapply(sites_core,length))
      site_nsamples<-data.frame(nsample=table(meta(newt1)$site))
      site_nsamples$ncore<-site_ncore_taxa[,1][match(site_nsamples[,1],rownames(site_ncore_taxa))]
      ggplot(site_nsamples,aes(x=nsample.Freq,y=ncore)) + geom_smooth(method = "lm") + geom_point(shape=21,size=5,fill="white") + theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

      site_ncore_taxa$site<-rownames(site_ncore_taxa)
      
### Strip Out Taxonomy 
      site_core_family_prop<-lapply(lapply(sites_core_taxonomy,function(x) table(x$Family)),function(x) x/sum(x))
      site_core_family_prop2<-lapply(site_core_family_prop,function(x) data.frame(family=names(x),prop=x))
      for(i in 1:length(site_core_family_prop2)){site_core_family_prop2[[i]]$site<-names(site_core_family_prop2)[i]}
      site_core_df<-do.call(rbind.data.frame, site_core_family_prop2)
      colnames(site_core_df)<-c("Family","Prop","Proportion","Site")
      corefamtab<-(table(site_core_df$Family)); corefams<-names(corefamtab)[corefamtab>1]
      site_core_df$Family_collapse<-ifelse(site_core_df$Family %in% corefams,site_core_df$Family,"Other")
      site_core_df$Family_collapse<-factor(site_core_df$Family_collapse,levels=c("Burkholderiaceae","Comamonadaceae","Oxalobacteraceae","Pseudoalteromonadaceae","Other"))

### Plot
      site_core_cols<-c(brewer.pal(5,"Paired"))
      site_core_plot1<-ggplot(site_core_df) + geom_bar(aes(x=Site,y=Proportion,fill=Family_collapse),stat="identity",color="grey40") + scale_fill_manual(values=site_core_cols) + theme_bw(base_size=15) + labs(fill="Family") + theme(legend.position = "top") #+ guides(fill=guide_legend(nrow=2))
      site_core_plot2<-site_core_plot1 + geom_text(data=site_ncore_taxa,aes(x=site,y=1,label=ntaxa),vjust=-0.2,cex=8) + guides(fill=guide_legend(nrow=3))
      site_core_plot2

      ggsave2('Core Visit 1 by Site.pdf',site_core_plot2,width=30,height=18,units="cm")

Disease Venn

### Plot Venn
  core_venn_cols<-brewer.pal(5,"Set2")[c(1:4)]
  core_vennplot1<-ggvenn(sites_core,fill_color = core_venn_cols,text_size = 8,set_name_size = 10)
  core_vennplot1

Combined Core Plot

core_grid<-plot_grid(site_core_plot2,core_vennplot1,labels="AUTO",ncol=2,label_size = 30)
ggsave2('Fig_6_Newt Core Microbiomes.pdf',core_grid,width=45,height=25,units="cm")

Genus Level Core Taxonomy

### Strip Out Taxonomy 
      site_core_genus<-lapply(lapply(sites_core_taxonomy,function(x) table(x$Genus)),function(x) x/sum(x))
      site_core_genus2<-lapply(site_core_genus,function(x) data.frame(genus=names(x),prop=x))
      for(i in 1:length(site_core_genus2)){site_core_genus2[[i]]$site<-names(site_core_genus2)[i]}
      site_core_genus_df<-do.call(rbind.data.frame, site_core_genus2)
      colnames(site_core_genus_df)<-c("Genus","Genus","Proportion","Site")
      write.csv(site_core_genus_df,'Site Core ASVs Genus.csv')

GENERAL LINEAR LATENT VARIABLE MODELS

Models

#Infected Sites
newt_infected<-prune_samples(sample_data(newt1)$site %in% c("R4","Rnew"),newt1)

## OTUs Top 50
  newt_infected_genus<-aggregate_top_taxa2(newt_infected, "Genus", top = 50)
  clr_scaled <-microbiome::transform(newt_infected_genus, transform = "clr")

#Extract
  ys <- data.frame(t(otu_table(clr_scaled)))
  names(ys) <-taxa_names(clr_scaled)
  
#Predictors
  Xs<-data.frame(sample_data(clr_scaled)) %>% dplyr::select(site,infection,svl)
 
## Model 
  fit_reduced_scaled <- gllvm(ys, Xs, 
                              num.lv = 2,
                              formula = ~ site * infection,
                              family = "gaussian",starting.val='zero')
  
    coefplot(fit_reduced_scaled)

    #plot(fit_reduced_scaled)

Plots of Effects

#Estimates 
  df<-coef(fit_reduced_scaled)
  est_df<-data.frame(df$Intercept)
  est_df2<-data.frame(df$Xcoef) 
  est_df3<-merge(est_df, est_df2, by = 0)
  
#Order genera
  row.names(est_df3)<-est_df3$Row.names
  est_df3<-est_df3[colnames(ys),]
  names(est_df3)[1]<- "Genus"
  names(est_df3)[2]<- "Intercept"
  

### COnfidence Intervals
  confint_df<-data.frame(confint(fit_reduced_scaled))
  
#Identify Rows with Main Effects
   one_comma<-sapply(rownames(confint_df),function(x) length(gregexpr(":", x, fixed = TRUE)[[1]]))==1
   
###Strip Out Individual Datasets
  ##R4
    r4_main<-grepl("^Intercept",rownames(confint_df))
    int_data_r4<-cbind(est_df3[,c("Genus","Intercept")],confint_df[r4_main,])
  ##Rnew
    rnew_main<-grepl("^Xcoef.siteRnew",rownames(confint_df)) & one_comma==TRUE
    int_data_rnew<-cbind(est_df3[,c("Genus","siteRnew")],confint_df[rnew_main,])
  ##Infection
    infection_main<-grepl("Xcoef.infection",rownames(confint_df))
    int_data_infection<-cbind(est_df3[,c("Genus","infection")],confint_df[infection_main,])
  #Rnew:Infection
     rnew_infect_main<-grepl("^Xcoef.siteRnew:infection",rownames(confint_df))
    int_data_rnew_infect<-cbind(est_df3[,c("Genus","siteRnew.infection")],confint_df[rnew_infect_main,])

  
#Extra Variables and Rename Columns 
  colnames(int_data_rnew)<-c("Genus","Estimate","l95","u95")
  colnames(int_data_infection)<-c("Genus","Estimate","l95","u95")
  colnames(int_data_rnew_infect)<-c("Genus","Estimate","l95","u95")
 
  int_data_rnew$trait<-"Site Rnew"
  int_data_infection$trait<-"Infection"
  int_data_rnew_infect$trait<-"Rnew:Infection"
    
  newt_mod_plotdata<-rbind(int_data_infection,int_data_rnew_infect)

#Order   
  newt_mod_plotdata$trait<-factor(newt_mod_plotdata$trait,levels=c("Infection","Rnew:Infection"))
  
  newt_mod_plotdata2<- newt_mod_plotdata %>% group_by(trait)  %>% arrange(Estimate,.by_group=T)
  newt_mod_plotdata2$Genus<-factor(newt_mod_plotdata2$Genus,levels=unique(newt_mod_plotdata2$Genus))
#Significance  
  newt_mod_plotdata2$Sig<- !data.table::between(0, newt_mod_plotdata2$l95, newt_mod_plotdata2$u95)
   
  sig_col<-brewer.pal(5,"Set1")[2]

### plot  
  newt_plot1<-ggplot(newt_mod_plotdata2,aes(x=Estimate,y=Genus)) + geom_errorbar(aes(y=Genus,xmin=l95,xmax=u95,color=Sig),alpha=0.6) + geom_point(size=5,shape=21,color="gray40",aes(fill=Sig),alpha=0.9)
  
  newt_plot2<- newt_plot1 + theme_bw(base_size = 15) + geom_vline(xintercept=0,linetype="dashed") + scale_color_manual(values=c("gray40",sig_col)) + scale_fill_manual(values=c("white",sig_col)) + guides(fill="none",color="none") + facet_wrap(.~trait,scales = "free_x")
  newt_plot2

Correlation Plots

#### Correlations
  cr1<-data.frame(getResidualCor(fit_reduced_scaled))#
  
  names(cr1)<-names(ys)
  names(cr1)<-abbreviate(names(cr1), minlength = 15)
  rownames(cr1)<-abbreviate(rownames(cr1), minlength = 15 )
  
  library(rstatix)
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:MASS':
## 
##     select
## The following object is masked from 'package:stats':
## 
##     filter
  #devtools::install_github("kassambara/ggcorrplot")
  library(ggcorrplot)
## 
## Attaching package: 'ggcorrplot'
## The following object is masked from 'package:rstatix':
## 
##     cor_pmat
  #install.packages("ggpubr")
  library(ggpubr)
## 
## Attaching package: 'ggpubr'
## The following object is masked from 'package:cowplot':
## 
##     get_legend
  cr2<-cor_pmat(cr1)
  
  corplot<-ggcorrplot(cr1, 
                      hc.order = TRUE,
                      outline.col = "white",
                      type = "full",
                      ggtheme = ggplot2::theme_minimal(base_size = 10),
                      tl.cex = 12,
                      p.mat = cr2,
                      sig.level = 0.05,
                      lab_size = 15,
                      #show.diag = F,
                      insig = "blank",
                      # colors = c("#6D9EC1", "white", "#E46726"))
                      colors = c("blue", "white", "red"))+
    theme(axis.text.x = element_text(angle = 90, vjust=0.5))+
    theme(plot.margin=unit(c(0.2,0.2,0.2,2),"cm"))
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the ggcorrplot package.
##   Please report the issue at <https://github.com/kassambara/ggcorrplot/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
  corplot

## Save Grid
  #newt_mod_combined1<-ggarrange(newt_plot2,corplot,nrow=2,widths=c(1,0.8),align="h",labels="AUTO",font.label=list(size=30))
  #newt_mod_combined1
  cowplot::ggsave2('Newt GLLVM Residual Correlations.pdf',corplot,width=25,height=25,units = "cm")

PREDICTED FUNCTIONAL ANALYSIS

Prepare for MicFunPred

Adapted from https://forum.qiime2.org/t/importing-dada2-and-phyloseq-objects-to-qiime-2/4683

# ## Trim Newt1
# newt1_filter<-prune_taxa(taxa_sums(newt1)>200,newt1)
# 
# #Make Fasta File From Sequences
# #remotes::install_github("alexpiper/seqateurs")
# library(seqateurs)
# ps_to_fasta(newt1_filter)
# 
# #Export Feature Table 
#   #library(biomformat);packageVersion("biomformat")
#    otu<-t(as(otu_table(newt1_filter),"matrix")) # 't' to transform if taxa_are_rows=FALSE
#   #Rename ASVs to Match FASTA above
#    rownames(otu)<-paste0("ASV_",seq(nrow(otu)))
#   #Fix Sample Names
#    true_samples_newt1filter<-colnames(otu)
#    colnames(otu)<-paste0("sample",seq(ncol(otu)))
#   #Write File 
#     # otu_biom<-make_biom(data=otu)
#     # write_biom(otu_biom,"Newt1_filter.biom")
#    write.table(otu,'newt1_filter.tsv',sep="\t",row.names=TRUE, col.names=TRUE, quote=FALSE)
# 
# ## Sample Data
# newt1_samplemeta<-as(sample_data(newt1_filter),"data.frame")
# write.csv(newt1_samplemeta,"Newt1_Filter_sample-metadata.csv")

BETA DIVERSITY OVER TIME

# Visit 1 and 2
  newt2<-prune_samples(sample_data(ps_clean_filter)$visit %in% seq(2) & !sample_data(ps_clean_filter)$site %in% c("R44","Rnew"),ps_clean_filter)

##CLR Transform - Round 1 & 2 sampling only
  newt2_clr <- microbiome::transform(newt2, "clr") 

##Extract Matrix and Sample Data
  newt2_clr_v<-vegan_otu(newt2_clr) #???
  newt2_clr_s<-as(sample_data(newt2_clr),"data.frame")

###Principal Components Analysis 
  newt2_pc<-prcomp(newt2_clr_v)
  summary(newt2_pc)$importance[,1:2]
##                             PC1      PC2
## Standard deviation     15.37126 14.99021
## Proportion of Variance  0.16088  0.15301
## Cumulative Proportion   0.16088  0.31389
## Metadata  
  newt2_pc_scores<-scores(newt2_pc)
  newt2_pc_scores_sub<-newt2_pc_scores[,1:2]
  newt2_pc_scores_sub<-cbind(newt2_pc_scores_sub,newt2_clr_s)

##Housekeeping + Label Setup    
  newt2_pc_scores_sub$site<-as.factor(newt2_pc_scores_sub$site)
  newt2_pc_scores_sub$visit<-as.factor(newt2_pc_scores_sub$visit)
  newt2_pc_scores_sub$site_visit<-with(newt2_pc_scores_sub,paste(site,visit))
  
### Plot  Data
  newt2plotdata<-gg_ordiplot(newt2_pc,groups=newt2_pc_scores_sub$site_visit,spider=T,plot=T)

  #newt2plotdata$plot + scale_color_brewer(palette = "Paired") + theme_bw() + geom_point(size=5)
  
  newt2plotdata2<-newt2plotdata$df_ord
  newt2plotdata3<-cbind(newt2plotdata2,newt2_clr_s)
  newt2plotdata3$Group<-factor(newt2plotdata3$Group)
  newt2plotdata3$infection_text<-ifelse(newt2plotdata3$infection==0,"No Visible Infection","Infected")
  
  gg_twovisits1<-ggplot() + geom_segment(data = newt2plotdata$df_spiders, aes(x = cntr.x, xend = x, y = cntr.y, yend = y, color = Group), 
    show.legend = FALSE) + geom_point(data=newt2plotdata3,aes(x=x,y=y,fill=Group,shape=infection_text),color="white",size=5)
  gg_twovisits2<- gg_twovisits1 + geom_path(data = newt2plotdata$df_ellipse, aes(x = x, y = y, color = Group), show.legend = FALSE)
  gg_twovisits3<- gg_twovisits2 + theme_bw(base_size=15)  + theme(legend.position = "right") + scale_shape_manual(values=c(24,21))
  gg_twovisits4<- gg_twovisits3 + guides(shape = guide_legend(override.aes = list(color = "black") ) ) + labs(x="PC1 (16.09%)",y="PC2 (15.3%)",fill="Site  & Visit",shape="Amphibiothecum \n Infection") + guides(fill = guide_legend(override.aes = list(shape=21))) 
  gg_twovisits5<-gg_twovisits4 + scale_color_brewer(palette="Paired") + scale_fill_brewer(palette = "Paired")
  gg_twovisits5

  ggsave2('Newt 2 Time Point PCa.pdf',gg_twovisits5,width=24,height=15,units="cm")

PERMANOVA

##Run PERMANOVA on CLR transformed data
visit_permanova<- adonis2(newt2_clr_v ~ site*visit + infection,data=newt2_clr_s,nperm=999,method="euclidean")
visit_permanova #
## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = newt2_clr_v ~ site * visit + infection, data = newt2_clr_s, method = "euclidean", nperm = 999)
##           Df SumOfSqs      R2      F Pr(>F)    
## Model      6    61508 0.37731 10.604  0.001 ***
## Residual 105   101509 0.62269                  
## Total    111   163017 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1